#define USE_PUPY_WEIGHTS true

#include <assert.h>
#include "nw_aligner_simple.cpp"

const char ALPHABET[]={'a', 'c', 't', 'g'};
#define ALPHABET_SIZE 4

struct MKResult {
    double ev;
    double sd;
};


MKResult run_monte_carlo(string seq1, string seq2, int iter_num) {
    assert (seq1.length() == seq2.length());
    int length = seq1.length();
    MICell* matrix = new MICell[(length+1)*(length+1)];
    
    int i, rand_score;
    double scores_sum = 0;
    double diff_sum = 0;
    double ev, sd;
    vector<int> scores(iter_num);
    
    // Monte-Carlo
    for (i=0; i<iter_num; i++) {
        random_shuffle(seq1.begin(), seq1.end());
        rand_score = _nw_score_restricted(seq1, seq2, matrix);
        scores_sum += rand_score;
        scores[i] = rand_score;
        //cout << "MK iter=" << i << "\t score=" << rand_score << "\n";
    }
    delete [] matrix;
    
    // calculating Z
    ev = scores_sum / iter_num;
    for (i=0; i<iter_num; i++) { diff_sum += (scores[i] - ev)*(scores[i] - ev); }
    sd = sqrt(diff_sum/iter_num);
    
    return {ev, sd};
}

extern "C"
MKResult run_monte_carlo_c(char *seq1_c, char *seq2_c, int iter_num) {
    string seq1 = (const char*) seq1_c;
    string seq2 = (const char*) seq2_c;

    return run_monte_carlo(seq1, seq2, iter_num);
}



int main( int argc, char *argv[] ) {
    int start_pos, length, prev_length, i, rand_score, score;
    double scores_sum;
    double ev, sd, diff_sum;
    int iter_num = 100000;
    vector<int> scores(iter_num);
    
    int min_length = 100;
    int max_length = 10000;
    float step_mult = 1.05;
    
    MICell* matrix = new MICell[(10000+1)*(10000+1)];

    string seq1 = "";
    string seq2 = "";

    for (length=min_length; length<=max_length; length*=step_mult) {
        
        prev_length = seq1.length();
        //for (i=0; i<ALPHABET_SIZE; i++) {
        //    seq1.append((length-prev_length)/ALPHABET_SIZE, ALPHABET[i]);
        //    seq2.append((length-prev_length)/ALPHABET_SIZE, ALPHABET[i]);
        //}

        seq1.append((length-prev_length)/10*4, ALPHABET[0]);
        for (i=1; i<ALPHABET_SIZE; i++) {
            seq1.append((length-prev_length)/10*2, ALPHABET[i]);
        }
        
        seq2.append((length-prev_length)/10*4, ALPHABET[0]);
        for (i=1; i<ALPHABET_SIZE; i++) {
            seq2.append((length-prev_length)/10*2, ALPHABET[i]);
        }

//        cout << seq1 << "\n";
//        cout << seq2 << "\n";

        //MKResult mk_res = run_monte_carlo(seq1, seq2, iter_num);
        //cout << mk_res.ev << '\t' << mk_res.sd << '\n';
        
        if (true) {    
 
            scores_sum = 0;
            
            // Monte-Carlo
            for (i=0; i<iter_num; i++) {
                random_shuffle(seq1.begin(), seq1.end());
                rand_score = _nw_score_restricted(seq1, seq2, matrix);
                scores_sum += rand_score;
                scores[i] = rand_score;
                //cout << "MK iter=" << i << "\t score=" << rand_score << "\n";
            }
            
            // calculating Z
            ev = scores_sum / iter_num;
            diff_sum = 0;
            for (i=0; i<iter_num; i++) { diff_sum += (scores[i] - ev)*(scores[i] - ev); }
            sd = sqrt(diff_sum/iter_num);
        
            sort(scores.begin(), scores.end());

            printf("%7i,%i", length, seq1.length());
            printf("\t%8.1f\t%8.1f", ev, sd);
            score = scores[0];
            printf("\t%i\t%.2f", score, (score - ev) / sd);
            score = scores[iter_num/10000*9500];
            printf("\t%i\t%.2f", score, (score - ev) / sd);
            score = scores[iter_num/10000*9900];
            printf("\t%i\t%.2f", score, (score - ev) / sd);
            score = scores[iter_num/10000*9990];
            printf("\t%i\t%.2f", score, (score - ev) / sd);
            score = scores[iter_num/10000*9999];
            printf("\t%i\t%.2f", score, (score - ev) / sd);
            score = scores[iter_num-1];
            printf("\t%i\t%.2f", score, (score - ev) / sd);
            printf("\n");
            
        };
    }

    return 0;

}
